
import cartopy.crs as crs
import cartopy.io.shapereader as shpreader
from cartopy.feature import ShapelyFeature, NaturalEarthFeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt
import numpy as np
import pygrib as pb
import matplotlib.colors as mcolors
import sys, glob
import matplotlib.ticker as mticker

datapath = glob.glob("./M-*grb2")
print(datapath)

['./M-A0064-21072000-000.grb2', './M-A0064-21072000-006.grb2', './M-A0064-21072000-012.grb2', './M-A0064-21072000-018.grb2', './M-A0064-21072000-024.grb2', './M-A0064-21072000-030.grb2', './M-A0064-21072000-036.grb2', './M-A0064-21072000-042.grb2', './M-A0064-21072000-048.grb2', './M-A0064-21072000-054.grb2', './M-A0064-21072000-060.grb2', './M-A0064-21072000-066.grb2', './M-A0064-21072000-072.grb2', './M-A0064-21072000-078.grb2', './M-A0064-21072000-084.grb2']


"""
利用cartopy套件繪製圖形，在read_grib範例中，可以得知變數的位置，這邊畫雨量跟10米風，初始場沒有雨量值
需要的步驟，
1.讀取grib資料，取出所需資料，包含2米溫度、10米u風、10米v風、經緯度
2.投影相關設定值(lon_0, lat_0, lat_1, lat_2)
  lon_0的key:LoVInDegrees
  lat_0的key:LaDInDegrees
  lat_1的key:Latin1InDegrees
  lat_2的key:Latin2InDegrees
3.繪製圖形
  這邊會先設定雨量值的colorbar
"""
cwb_data=['None','#9BFFFF','#00CFFF','#0198FF','#0165FF','#309901','#32FF00','#F8FF00','#FFCB00','#FF9A00',\
          '#FA0300','#CC0003','#A00000','#98009A','#C304CC','#F805F3','#FECBFF']
clevs = [0,1,2,6,10,15,20,30,40,50,70,90,110,130,150,200,300,400]
cmap = mcolors.ListedColormap(cwb_data,'preecipitation')
norm = mcolors.BoundaryNorm(clevs, cmap.N)



cnt = 0 # 只是為了抓取1測的固定值用 
for grb in range(len(datapath)-1):
    if cnt == 0:
        forget = pb.open(datapath[grb]).select()[:][0]
        analDate = forget.analDate.strftime("%y%m%d%H")
        lon_0, lat_0 = forget["LoVInDegrees"], forget["LaDInDegrees"]
        lat_1, lat_2 = forget["Latin1InDegrees"], forget["Latin2InDegrees"]
        lats, lons = forget.latlons()
        proj = crs.LambertConformal(central_longitude=lon_0, central_latitude=lat_0,standard_parallels=(lat_1,lat_2))
        coast = NaturalEarthFeature(category='physical', scale='10m',
                            facecolor='none', name='coastline')
    #gispath ="/mnt/c/cwbconda/shapefile"
    basicvar1 = pb.open(datapath[grb]).select()[:]
    basicvar2 = pb.open(datapath[grb+1]).select()[:]
    fcst = str(basicvar2[0].forecastTime)
    valDate = basicvar2[0].validDate.strftime("%y%m%d%H")
    accrain = basicvar2[61]["values"] - basicvar1[61]["values"]
    u10m_data = basicvar2[66]["values"]
    v10m_data = basicvar2[67]["values"]
    fig = plt.figure(figsize=(20,16))
    axs = plt.axes(projection=proj)
    axs.add_feature(coast, edgecolor='gray')
    ctf = axs.contourf(lons,lats, accrain, clevs,cmap=cmap, norm=norm, transform=crs.PlateCarree())
    cbar = plt.colorbar(ctf,ticks=clevs[1:len(clevs)-1],orientation='horizontal',drawedges=True,pad=0.05)
    cbar.ax.tick_params(labelsize=13)
    cbar.set_label("mm",size=14)
    cbar.dividers.set_color('white')
    cbar.dividers.set_linewidth(5)
    cbar.outline.set_color('white')
    #cbar.outline.set_linewidth(2)
    gap = 25
    axs.barbs(lons[::gap,::gap],lats[::gap,::gap],u10m_data[::gap,::gap],v10m_data[::gap,::gap],length=6,\
               color="midnightblue",lw=1.5,transform=crs.PlateCarree())
    gl = axs.gridlines(draw_labels=True,x_inline=False, y_inline=False)
    gl.xlabels_top= False
    gl.ylabels_right = False
    #gl.xlocator = mticker.FixedLocator([118,120,122,124,126,128])
    #gl.ylocator = mticker.FixedLocator([19,21,23,25,27])
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    axs.set_extent([lons[0,0], lons[0,-1], lats[0,0], lats[-1,0]])
    axs.set_title("Initial Time: {}".format(analDate), loc="left",fontsize=16)
    axs.set_title("Forecast Time:{}".format(fcst),loc="center",fontsize=16)
    axs.set_title("Valid Time:{}".format(valDate),loc="right",fontsize=16)
    plt.suptitle("In-Fa Typhoon",y=0.91,fontsize=20)
    plt.show()
    plt.clf()

/home/louis/anaconda3/envs/summer/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:307: UserWarning: The .xlabels_top attribute is deprecated. Please use .top_labels to toggle visibility instead.
  warnings.warn('The .xlabels_top attribute is deprecated. Please '
/home/louis/anaconda3/envs/summer/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:343: UserWarning: The .ylabels_right attribute is deprecated. Please use .right_labels to toggle visibility instead.
  warnings.warn('The .ylabels_right attribute is deprecated. Please '

/home/louis/anaconda3/envs/summer/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:307: UserWarning: The .xlabels_top attribute is deprecated. Please use .top_labels to toggle visibility instead.
  warnings.warn('The .xlabels_top attribute is deprecated. Please '
/home/louis/anaconda3/envs/summer/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:343: UserWarning: The .ylabels_right attribute is deprecated. Please use .right_labels to toggle visibility instead.
  warnings.warn('The .ylabels_right attribute is deprecated. Please '

<Figure size 432x288 with 0 Axes>

/home/louis/anaconda3/envs/summer/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:307: UserWarning: The .xlabels_top attribute is deprecated. Please use .top_labels to toggle visibility instead.
  warnings.warn('The .xlabels_top attribute is deprecated. Please '
/home/louis/anaconda3/envs/summer/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:343: UserWarning: The .ylabels_right attribute is deprecated. Please use .right_labels to toggle visibility instead.
  warnings.warn('The .ylabels_right attribute is deprecated. Please '

<Figure size 432x288 with 0 Axes>

/home/louis/anaconda3/envs/summer/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:307: UserWarning: The .xlabels_top attribute is deprecated. Please use .top_labels to toggle visibility instead.
  warnings.warn('The .xlabels_top attribute is deprecated. Please '
/home/louis/anaconda3/envs/summer/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:343: UserWarning: The .ylabels_right attribute is deprecated. Please use .right_labels to toggle visibility instead.
  warnings.warn('The .ylabels_right attribute is deprecated. Please '

<Figure size 432x288 with 0 Axes>

/home/louis/anaconda3/envs/summer/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:307: UserWarning: The .xlabels_top attribute is deprecated. Please use .top_labels to toggle visibility instead.
  warnings.warn('The .xlabels_top attribute is deprecated. Please '
/home/louis/anaconda3/envs/summer/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:343: UserWarning: The .ylabels_right attribute is deprecated. Please use .right_labels to toggle visibility instead.
  warnings.warn('The .ylabels_right attribute is deprecated. Please '

<Figure size 432x288 with 0 Axes>

/home/louis/anaconda3/envs/summer/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:307: UserWarning: The .xlabels_top attribute is deprecated. Please use .top_labels to toggle visibility instead.
  warnings.warn('The .xlabels_top attribute is deprecated. Please '
/home/louis/anaconda3/envs/summer/lib/python3.7/site-packages/cartopy/mpl/gridliner.py:343: UserWarning: The .ylabels_right attribute is deprecated. Please use .right_labels to toggle visibility instead.
  warnings.warn('The .ylabels_right attribute is deprecated. Please '